function fig4_10

global A

clf
% get(gcf)
set(gcf,'Position', [25 1115 658 230])
co = [0 0 1;
    0 0.5 0;
    1 0 0;
    0 0.75 0.75;
    0.75 0 0.75;
    0.75 0.75 0;
    0.25 0.25 0.25];
set(groot,'defaultAxesColorOrder',co)
subaxis(1,1,1,1,'MT',0.01,'MB',0.2,'MR',0.01,'ML',0.05,'P',0.04)
hold on

tmax=800;

m1=1; m2=1; k1=1;

k2=k1;
%k2=(1.1)^2*k1;

w2=sqrt(k2/k1)
k12=1; c12=1;
ep=0.01;

if k2==k1
    say=['$\omega_2 = \omega_1$'];
else
    say=['$\omega_2 = 1.1\omega_1$'];
end
text(600,0.75,say,'Interpreter','latex','FontSize',16)

km1=k1/m1; km2=k2/m2;
em1=ep/m1; em2=ep/m2;

A=zeros(4,4);
A(1,2)=1;
A(2,1)=-km1-em1*k12; A(2,2)=-em1*c12; A(2,3)=em1*k12; A(2,4)=em1*c12;
A(3,4)=1;
A(4,1)=em2*k12; A(4,2)=em2*c12; A(4,3)=-km2-em2*k12; A(4,4)=-em2*c12;

y0=zeros(4,1);
y0(1)=1;
y0(3)=-0.5;

%  calculate solution using a MATLAB routine
[t,y] = ode45(@rhs,[0 tmax],y0);
%[t,y] = ode23s(@rhs,[0 tmax],y0);  titleX='ode23s';

hold on

plot(t,y(:,3),'r','LineWidth',1)
plot(t,y(:,1),'b','LineWidth',0.5)
grid on
box on

set(gca,'ytick',[-1 0 1])
set(gca,'YTickLabel',{'$-1$';'0';'1'})

ax = gca;
ax.XAxis.TickLabelInterpreter = 'latex';
ax.YAxis.TickLabelInterpreter = 'latex';
set(gca,'FontSize',13)

ylabel('$y$\,-\,axis','Interpreter','latex','FontSize',14)
xlabel('$t$\,-\,axis','Interpreter','latex','FontSize',14)

% if k2==k1
%     exportgraphics(gcf,'/Users/mark/Desktop/oscil1.pdf')
% else
%     exportgraphics(gcf,'/Users/mark/Desktop/oscil2.pdf')
% end

%  define f(t,y)
function dy=rhs(t,y)
global A
dy=A*y;








































